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Abstract. - In order to clarify how the statistical properties of earthquakes depend on 
the constitutive law characterizing the stick-slip dynamics, we make an extensive numerical 
simulation of the one-dimensional spring-block model with a rate- and state-dependent friction 
law. Both the magnitude distribution and the recurrence-time distribution are studied with 
varying the constitutive parameters characterizing the model. While a continuous spectrum 
of seismic events from smaller to larger magnitudes is obtained, earthquakes described by this 
model turn out to possess pronounced "characteristic" features. 



Recent studies have revealed that an earthquake can be regarded as a stick-slip frictional 
instability which a natural fault driven by steady motions of tectonic plates exhibits. Hence, an 
earthquake occurrence is governed by the physical law of rock friction [1,2]. Unfortunately, our 
present understanding of physics of friction is still poor. We do not have precise knowledge of 
the constitutive law characterizing the stick-slip dynamics of earthquake faults. The difficulty 
lies partly in the fact that a complete microscopic theory of friction is still not available, but 
also in the fact that the length and time scales relevant to earthquakes are so large that the 
applicability of laboratory experiments on rock friction is not necessarily clear. 

Detailed characteristics of the friction force are specified by the constitutive relation [1-4] . 
One fundamental question in earthquake studies might be how the properties of earthquakes 
depend on the constitutive law and other material parameters characterizing earthquake faults. 
To answer this question and to get deeper insight into the physical mechanism governing the 
stick-slip process of earthquakes, proper modeling of an earthquake is essential. Indeed, earth- 
quake models of various levels of simplifications have been proposed in geophysics and statis- 
tical physics, and their statical properties have been extensively studied mainly by means of 
numerical computer simulations. In model simulations, one can easily control the constitutive 
parameters, which is almost impossible for natural faults. 

One of the standard model used in statistical studies of earthquakes is the so-called spring- 
block model originally proposed by Burridge and Knopoff (BK model) [5]. The BK model, 
combined with several types of constitutive relations, have been extensively studied by numer- 
ical simulations [6-13]. In order for the model to exhibit a dynamical instability corresponding 
to an earthquake, it is essential that the friction force exhibits a frictional weakening property, 
i.e., the friction should become weaker as the block slides. One of the simplest form of the 
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friction force widely used is a velocity- weakening friction force [7, 8] . Here, the friction force 
is assumed to be a single- valued function of the velocity, getting smaller as the velocity in- 
creases. The other form employed in the previous analyses might be a slip- weakening friction 
force where the friction force is assumed to be a single- valued function of the cumulative fault 
slip, getting smaller as the slip distance increases [9, 10, 14]. 

From laboratory experiments of rock friction, however, there is an indication that the 
real constitutive relation might be more complex, neither purely velocity-weakening nor slip- 
weakening. Experiments suggest that the friction force depends on some "hidden" variable, 
possibly representing the state of the rock interface. Some time ago, Dieterich and Ruina 
proposed an empirical form of the constitutive law, a rate- and state-dependent friction law, 
by phcnomenologically introducing the time-dependent "state variable" and its time-evolution 
equation [3,4]. Though introduced empirically, this rate- and state-dependent constitutive 
law is devised so as to reproduce certain noticeable features of rock experiments mentioned 
above. Tse and Rice employed this rate- and state-dependent constitutive relation in their 
numerical simulations of earthquakes [15]. These authors studied the stick-slip motion of 
the two-dimensional strike-slip fault within an elastic continuum theory, assuming that the 
fault motion is rigid along strike. It was then observed that large events repeated quasi- 
periodically. Since then, similar rate- and state-dependent constitutive laws have widely been 
used in numerical simulations [16-21]. Somewhat different type of slip- and state-dependent 
constitutive law was also used [22]. 

Cao and Aki performed a numerical simulation of earthquakes by combining the one- 
dimensional BK model with the rate- and state-dependent friction law in which various con- 
stitutive parameters were set nonuniform over blocks [23]. In the present letter, we wish to 
extend an earlier calculation by Cao and Aki to study the statistical properties of the ID BK 
model combined with the rate- and state-dependent constitutive law with uniform constitutive 
parameters. Such a calculation might be useful due to the following two reasons: First, due to 
the simplicity of the ID BK model, it is now possible to generate large number of events for 
large enough system to obtain various statistical properties of sufficient accuracy. Second, by 
comparing the obtained results with the previous ones for the ID BK model with the purely 
velocity-weakening or slip-weakening constitutive law with uniform constitutive parameters, 
it is possible to clarify the dependence of earthquake properties on the underlying constitutive 
law. 

In the one-dimensional BK model, an earthquake fault is simulated by an array of blocks, 
each of which is connected to the neighboring blocks via the elastic springs of the spring 
constant k c , and to the moving plate via the springs of the spring constant k p . All blocks are 
subject to the friction force <f), the source of the nonlinearity in the model. 

The equation of motion for the i-th block can be written as 

mu'i = k v {v't' - u[) + k c {u' i+1 - 2u' t + u'^) - <j>(u'i,%), (1) 

where t' is the time, u[ the displacement of the i-th block, m the mass of the block, and v' 
is the loading rate representing the speed of the plate. The form of the friction force <j> is 
specified by the constitutive relation, which, as mentioned above, is a vitally important, yet 
largely ambiguous part in the proper description of earthquakes. Here, as the form of the 
friction force, we assume a rate- and state-dependent friction force, 

= {c' + a' log(l + % + b' log (2) 
v* L 

where v[ — u\ is the velocity of the i-th block, 9'^t) is the time-dependent state variable 
(with the dimension of time) representing the "state" of the slip interface, v* is a reference 
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velocity, N is an effective normal load, £ is a characteristic slip distance which is a measure 
of the distance of sliding necessary for the surface to evolve to a new state, while a' , b' and 
d are numerical constants describing the rate- and state-dependent friction law. The first 
term (c-term) is a constant taking a value around -|, which dominates the total friction in 
magnitude [1]. The second term (a-term) is a velocity-strengthening direct term describing 
the part of the friction which follows the velocity-change immediately. Note that we put a 
factor unity in the logarithm of this second term, which enables one to describe the system at a 
complete halt, whereas, without this term, the system cannot stop because of the logarithmic 
anomaly occurring at v' = 0. The third part (&-term) is an indirect term dependent of the 
state variable and follows the velocity-change via the time dependence of the state variable. 
Laboratory experiments suggest that the a- and 6-terms are smaller than the c-term by one 
or two orders of magnitudes [1,2], yet they play an essential role in stick-slip dynamics. 

There are several proposals for the form of the time-evolution equation of the state variable. 
Here, we assume the so-called slowness law given by 

We take the length unit to be the characteristic slip distance C, and the time unit to be 
the rise time of an earthquake w" 1 = (m/kp) 1 / 2 , and set v* = Clo. The coupled equations 
of motion are then made dimensionless by introducing the dimensionless variables, t — cut' , 
Ui = u'JC, Vi = v'Jv*, 9i = B'y/C, v = v'/v*, a = a'N/{k p C), b = b'N/{k p C), and 
c = c'M/(k p C), 

d 2 Ui 

= (vt - Ui) + l 2 {u i+ i - 2ui + Ui-i) - (c + alog(l + Vi) + blogOi) (4) 
= 1-vA (5) 



dt 2 
dfk 
dt 



where I = (kc/kp) 1 / 2 is the dimensionless stiffness parameter. 

With natural faults in mind, we give here rough estimates of various model parameters. 
Via the rise time of large earthquakes, the time unit uu^ 1 may be estimated to be a few 
seconds. Via the rupture-propagation speed, the block size d may be estimated to be a 
few kilometers [12]. The estimate of the characteristic slip distance C remains to be largely 
ambiguous: Here, we use an estimate by Scholz [1] and by Tse and Rice [15], C being of order 
a few mm or cm. Since the loading rate associated with the plate motion is typically a few 
cm/year, the dimensionless loading rate v = v'/(Cu) is of order v ~ 10~ 8 . The dimensionless 
quantity k p C/N may be written in terms of the normal stress a n , the density of the crust 
p and other parameters defined above as pduj 2 C/ a ni which may be estimated of order 10~ 4 . 
The dimensionless parameter c is then estimated to be of order 10 3 ~ 10 4 . As mentioned, the 
parameters a and b are one or two orders of magnitude smaller than c. In the following, we 
set the parameter values of the model around the typical values estimated above. 

The coupled equations of motion (2) are solved numerically by using the Runge-Kutta 
method of the fourth order. Total number of 2 x 10 5 events are generated in each run, which 
are used to perform various averagings, while initial 10 4 events are discarded as transients. 
The total number of blocks N are taken to be N = 800, with the open boundary condition. 
In some cases, larger systems with N = 1600 are simulated to check the size dependence. In 
the present simulation, we fix the loading rate to v = 10~ 8 . The clastic parameter I is set to 
I = 3, which is the value extensively studied by Mori and Kawamura for the ID BK model 
with the Carlson-Langer velocity- weakening friction law [11]. 
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Fig. 1 - The magnitude distribution of seismic events. In the upper figure, the parameter b is varied 
in the range 20 < b < 50 with fixing c = 10 3 and a = 0. In the middle figure, the parameter a is 
varied in the range < a < 30 with fixing c = 10 3 and b = 50. In the lower figure, the parameter c is 
varied in the range 10 2 < c < 10 4 with fixing a = and b — 20. The other parameters are v = 10 _s , 
I — 3 and N = 800 for all cases. The magnitude distribution depends solely on 6, not on a and c. 
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In Figs.l(a)-(c), we show the magnitude distribution R(fi) of earthquakes calculated from 
our numerical simulations. The magnitude of an event, /i, is defined by jj, = log 10 (^ i Au,), 
where Am is the total displacement of the i-th block during a given event and the sum is taken 
over all blocks involved in the event [7]. In all cases studied, the calculated i?(/x) exhibits a 
continuous nontrivial distribution from smaller to larger magnitudes. 

In the upper panel of Fig.l, the parameter b is varied in the range 20 < b < 50 where the 
other constitutive parameters are fixed to c = 10 3 and a = 0. As can be seen from the figure, 
the magnitude distribution depends on the 6-valuc considerably. In the case of b = 20, the 
distribution exhibits a pronounced peak at a large magnitude, with a pronounced feature of 
a characteristic earthquake. In the case of b = 30 ~ 50, the distribution becomes flatter. In 
either case, the distribution is not a power-law, largely deviating from the Gutcnberg-Richtcr 
(GR) law. We also made a preliminary study of the smaller b region (b < 20) and of the larger 
b region (b > 50). In the case of smaller 6, an event tends to occur as a creep- like small event 
rather than a stick-slip, whereas, in the case of larger b, more weight is given to system-size 
large events in i?(/z) and the data suffer from significant finite-size effects. 

In the middle panel of Fig.l, the parameter a is varied in the range < a < 30 with 
c = 10 3 and b = 50. The magnitude distribution hardly depends on the a-value, indicating 
that the relavant parameter here is b itself, rather than a — b as widely believed [1,2]. In the 
lower panel of Fig.l, the parameter c is varied in the range 10 2 < c < 10 4 with a = and 
b = 20. Again, the magnitude distribution hardly depends on the c-value. Hence, it turns out 
that the magnitude distribution is sensitive only to the 6-value. 

We also calculate the local recurrence-time distribution of large events in order to see how 
large events repeat in time: In defining the recurrence time locally, the subsequent large event 
is counted when a large event occurs with its epicenter in the region within 10 blocks from 
the epicenter of the previous large event. In Figs. 2, we show the distribution of the local 
recurrence time T of large earthquakes whose magnitude is greater than fj, > fi c = 4 for the 
case of b = 20 (upper panel) and b = 40 (lower panel). The other constitutive parameters are 
fixed to c = 10 3 and a = 0. 

As can be seen from the upper panel of Fig. 2, the distribution in the case of b = 20 shows 
a sharp peak at Tv ~ 500, indicating a pronounced periodic feature of earthquake recurrence. 
Smaller peaks are also discernible at multiples of an elementary period. The inset represents 
a magnified view of the sharp peak at Tv ~ 500, where the data both for N = 800 and 1600 
are shown. The peak tends to get lower with N . In order to evaluate the possible boundary 
effect, we also calculate the recurrence-time distribution only for the subset of events which 
stay in the interior of the system, i.e., the events whose rupture does not reach the boundary, 
and the results are shown in the inset for both A?" = 800 and 1600. The data for these "interior 
events" no longer show appreciable size effect suggesting that this sharp peak might persist 
even in an infinite-size limit. 

Somewhat different behavior is observed in the case of b = 40. As can be seen from the 
lower panel of Fig. 2, the distribution exhibits two independent peaks at Tv ~ 400 and at 
Tv ~ 1000. Essentially the same behavior is also observed for b — 30 and 50. The inset 
represents a magnified view of the second peak at Tv ~ 1000. While the peak tends to get 
lower with increasing N, it is still clearly visible even at N = 1600. The first peak located 
at Tv ~ 400, by contrast, does not show appreciable size effect. The observed double-peak 
structure means that events tend to occur with double periods, although the origin of such a 
doubly-periodic behavior is not clear at the present stage. 

In summary, we studied the statistical properties of the ID BK model combined with 
a rate- and state-dependent constitutive law. We found the followings: (i) A continuous 
size distribution, spanning from smaller to larger magnitudes, has been obtained. Somewhat 
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Fig. 2 - The local recurrence-time distribution of large events of /i > fi c — 4 for the case of b = 20 
(upper panel) and of b = 40 (lower panel). The other parameters are c = 10 3 , a = 0, v — 10 -8 and 
/ = 3. The system size is N = 800 or 1600. The insets represent magnified views of the peak region. 
The data represented as "interior" are the ones taken for the type of large events whose rupture does 
not reach the boundary of the system. The observed distinct peak structure suggests near periodic 
occurrence of events. 



unexpectedly, however, the power-law distribution like the GR law is not realized in this 
model, in spite of the fact that the rate- and state-dependent friction law employed here is 
the standard one widely used in recent earthquake studies, (ii) Among several parameters 
characterizing the rate- and state-dependent constitutive law, the statistical properties of 
earthquakes depend most sensitively on the parameter b, which describes the extent of the 
state-dependent frictional instability, (iii) The recurrence-time distribution indicates the near 
periodic recurrence of large events. For the case of stronger frictional instability, doubly- 
periodic recurrence characterized by the two independent periods becomes eminent. 

These observations suggest that earthquakes described by the the present model possesses 
pronounced "characteristic" features. The characteristic tendency of our present model is even 
more enhanced than the one observed for the BK model with the Carlson-Langer velocity- 
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weakening constitutive law. Such pronounced characteristic features are in apparent contrast 
to the observations for natural faults, since the analysis of real seismic catalog is known to 
lead to more "critical" behaviors, e.g., the GR law for the magnitude distribution: See, e.g., 
Ref. [12]. 

Identifying the cause of this deviation would be of particular importance. First, we note 
that the rate- and state-dependent constitutive law was derived based on the data of rock 
friction in laboratory measurements. One has to be careful that the relevant length scale 
in laboratory measurements differs from the one at natural faults by orders of magnitude. 
There is no obvious guarantee that the rate- and state-dependent constitutive law derived 
from laberatory experiments is also valid in describing natural faults. Other possible reason 
of the deviation may be an intrinsic discreteness of the present model. According to Refs. 
[18, 19], however, the discreteness tends to enhance the criticality, which is just opposite to 
what is required to account for the observed deviation. We also need to recognize that real 
earthquake catalog is usually taken not for a single fault, but over many faults. There has been 
a suspicion that, even if the property of a single uniform fault is characteristic, the property 
obtained after averaging over many faults, each of which has different characteristics, might 
become apparently critical. If this is really the case, the observations for natural faults is not 
necessarily inconsistent with the observations for the present model, since the BK model deals 
with the property of a single uniform fault. In any case, further studies are required to settle 
the issue. 
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